Big Bang Nucleosynthesis with Independent Neutrino Distribution Functions 
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We have performed new Big Bang Nucleosynthesis calculations which employ arbitrarily-specified, 
time-dependent neutrino and antineutrino distribution functions for each of up to four neutrino 
flavors. We self-consistently couple these distributions to the thermodynamics, the expansion rate 
and scale factor-time/temperature relationship, as well as to all relevant weak, electromagnetic, 
and strong nuclear reaction processes in the early universe. With this approach, we can treat 
any scenario in which neutrino or antineutrino spectral distortion might arise. These scenarios 
might include, for example, decaying particles, active-sterile neutrino oscillations, and active-active 
neutrino oscillations in the presence of significant lepton numbers. Our calculations allow lepton 
numbers and sterile neutrinos to be constrained with observationally-determined primordial helium 
and deuterium abundances. We have modified a standard BBN code to perform these calculations 
and have made it available to the community. 

PACS numbers: 14.60.Pq; 14.60.St; 26.35,+c; 95.30.-k 



I. INTRODUCTION 



There is a new paradigm in Big Bang Nucleosyn- 
thesis (BBN) studies which promises enhanced probes 
of the early universe and a window into new physics. 
In the past, BBN predictions have been used to place 
constraints on the baryon number at three minutes af- 
ter the Big Bang. This was done by comparing the 
observationally-inferred primordial light element abun- 
dances to abundances predicted by BBN calculations 
over a wide range of baryon-to-photon ratio values. With 
the high precision results of the Wilkinson Microwave 
Anisotropy Probe (WMAP), however, the baryon-to- 
photon ratio, rj, is now independently determined - at 
300,000 years after the Big Bang - from observations of 
the cosmic microwave background (CMB) relative acous- 
tic peak amplitudes 0, H H- Currently, the WMAP 
Three Year Mean value for the baryon-to-photon ra- 
tio is 77 = (6.11 ± .22) x 10~ 10 . Future missions (e.g., 
Planck[|) promise considerably higher precision deter- 
minations Of 1]. 

Since the baryon-to-photon ratio is known indepen- 
dently, and to excellent precision albeit at much later 
times, BBN calculations can now be used to probe or 
constrain new physics or heretofore poorly determined 
parameters. For example, we can use BBN predictions 
to constrain not only the lepton numbers but also the 
physics behind these lepton numbers . The existence of a 
nonzero electron lepton number follows from charge neu- 
trality and the observed proton content of the universe. 
The contributions of neutrinos and antineutrinos to the 
electron, muon, and tau (e,/i, r) lepton numbers are not 
known, since we do not directly observe these relic parti- 
cles. The neutrino contribution to the lepton number for 
a given flavor, a = e, fi, r, is defined analogously to the 



baryon-to-photon ratio, rj = (rib — n l)/ n "ji as 



J1i/~ - 77^„ 



(1) 



where n 7 = (2£(3)/7r 2 )T^ is the proper photon number 
density at temperature T 7 , and n Va and n Pa are the neu- 
trino and antineutrino number densities. Observational 
bounds on the lepton numbers 0, i, 0, 0, ©, M, El El, El 
remain large compared to the values of these that could 
significantly affect BBN when there is new leptonic sector 
physics (e.g., sterile neutrinos) Q . 

The neutrino lepton numbers influence BBN and the 
resulting primordial element abundances in a number of 
wavs[14j. The energy density in the neutrino sector con- 
tributes to the total energy density of the universe which 
determines the expansion rate. The expansion rate is 
crucial to the outcome of BBN because it determines the 
weak freeze-out temperature which in turn effectively sets 
the neutron-to-proton ratio and, therefore, the primor- 
dial abundances of 4 He and the other light elements. 

Not only is the total number of neutrinos important 
to the outcome of BBN, but the neutrino distribution 
functions are key components of the phase space integrals 
in the weak reaction rates in BBN. The weak reactions 
of greatest interest are those that inter-convert neutrons 
and protons: 

v e + n^p + e~, (2) 



(3) 



-p + e +v e . 



(4) 



Since the rates for the weak reactions are strongly energy 
dependent, the energy distributions of the neutrinos and 
antineutrinos can figure prominently in both the forward 
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and reverse rates in the processes in Eqs. ([2]), ([3]), and 
(U|) . In standard BBN scenarios the neutrino distribution 
functions are assumed to be thermally-shaped Fermi- 
Dirac distributions. However, it is possible that non- 
thermal neutrino distribution functions arise after the 
neutrinos decouple from the background plasma around 
T w 3 MeV and during times crucial to BBN. 

There are many possible mechanisms that could alter 
the neutrino spectra. Altered neutrino energy spectra, 
in turn, could change the resulting primordial element 
abundances from what one would expect given a partic- 
ular lepton number. Neutrino energy spectrum-altering 
scenarios include, but are not limited to, active-active 



neutrino oscillations (E 



9l 1 1 011 , active-sterile neutrino 
oscillations [E [E 16 AYa, 1x81 ] , or particle decay into 
the neutrino sea[l9|. Moreover, active-sterile neutrino 
flavor mixing and other mechanisms for creating ster- 
ile neutrino dark matter before neutrino de coupling are 
a focus of current research[U H3, 113, HE H3, HE 
HI HI HE HE HI HI, as is the constraint of these sce- 
narios via x-ray observations and larg e-scale structure 
considerations [HEHIMSSSIlEliEliii- Though 
these models may not directly affect BBN through the 
spectral distortion of v e and v e energy distribution func- 
tions discussed here, they nevertheless may affect the 
overall values of lepton number, entropy, and energy den- 
sity which are relevant to BBN. In the end, the existence 
of sterile neutrino states changes the meaning and util- 
ity of lepton number 42, 43]. To use BBN predictions to 
probe or constrain any such scenario requires an approach 
that self-consistently includes neutrino and antincutrino 
energy spectra of arbitrary shape. 



We have performed detailed calculations of primor- 
dial nucleosynthesis in which we include neutrino and 
antineutrino spectral distortion. Our results are surpris- 
ing. We find that even modest distortions of the neu- 
trino and/or antineutrino spectral shapes from Fermi- 
Dirac black body forms can result in significant modi- 
fication of the net neutron-proton interconversion rates 
and, hence, alteration of the light element abundances. 

To study the effects of neutrino spectral distortion, 
we have modified the original Kawano/ Wagoner BBN 
code described in Ref. [44| to calculate the primordial 
element abundances self-consistently with arbitrarily- 
specified non-thermal and/or time-dependent neutrino 
distribution functions. This paper is structured as fol- 
lows: Section II describes the calculation of weak charge- 
changing reaction rates in the early universe and our 
prescription for employing non-thermal neutrino and an- 
tineutrino energy distribution functions; Section III dis- 
cusses our new BBN code; Section IV will present exam- 
ple results for non-thermal neutrino distribution func- 
tions resulting from various physical scenarios; and Sec- 
tion V gives conclusions. 




FIG. 1: Example neutrino occupation probabilities. The up- 
per dark (black) curve is the standard Fermi-Dirac thermally- 
distrubuted neutrino occupation probability and the lower 
light (red) curve is an example non-thermal neutrino occupa- 
tion probability which can result from active-sterile neutrino 
transformation. 



II. BBN AND THE WEAK REACTION RATES 

At early times and high temperatures, t ~ 1 sec and 
T > 1 MeV, the primordial element abundances are 
given by nuclear statistical equilibrium (NSE). In NSE 
the rates for the processes which create a particular nu- 
cleus are equal to the rates that destroy it, so that the 
abundance for each element is given by the Saha equa- 
tion. 

As the universe expands and cools, reaction rates slow 
down to the point where they will not be fast enough to 
maintain NSE and the neutron and proton abundances, 
and subsequently the abundances of 4 He and the other 
light nuclei, "freeze-out". For example, the 4 He abun- 
dance falls below its equilibrium NSE track at T m 0.6 
MeV, essentially as a consequence of the small NSE deu- 
terium abundance. BBN can be looked at crudely as 
a series of freeze-outs from NSE, but with considerable 
post-equilibrium nuclear processing. 

Because the entropy per baryon is high, alpha particles 
form copiously during BBN. Nearly all the neutrons in 
the universe at the epoch where a's form end up in alpha 
particles. 

A key factor in the outcome of BBN is the value of the 
neutron-to-proton ratio. Like the nuclear abundances in 
NSE, at high enough temperatures (T > 3 MeV) the weak 
neutron-proton inter-conversion rates are fast enough to 
maintain chemical equilibrium and the neutron-to-proton 
ratio can be determined from a Saha equation when the 
neutrinos have thermally-shaped distribution functions 
(as we will describe later). 

For general conditions the neutron-to-proton ratio is 
determined by the weak reaction processes shown in Eqs. 
([2HU). The rates for these weak reactions are given in 
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Eqs. (|5l fT0|) below. The forward rate for the reaction in 
Eq. ^ is given by X Uen , Eq. (jHJ), and the corresponding- 
reverse rate is given by A e - p , Eq. ([5|). Likewise, the for- 
ward and reverse rates for the process in Eq. ([3]) are \p eP 



and X e + n respectively. Eq. ([9]) gives the rate for free neu- 
tron decay denoted by A n _d 0C ay, while the reverse three- 
body reaction rate is denoted by A r p+ p. g iven in Eg . (jlOjl . 
These rates are detailed belowjl, |H, |M H3, |48L|49|: 



In 2 



(ft)(m e c 2 )' 



J F[Z,E u + Q n p]El{E v + Q np )[{E u + Q np ) 2 - m e c 2 j [S e -] [1 - S Ue ] dE v , 



(5) 



In 2 



(ft)(m e c 2 f 



/ E v (E v — Q np ) ( (E v — Q % 

JQ np +m c c 2 V 



V2 , 



[Sp e ] [1 - S e+ ] dE v 



(6) 



In 2 



(ft)(m e c 2 T 



% p +m B c- 



E*u {E v — Qnp) (y{E v — Qnp) 2 



1/2 



[S e +][l-S PB ]dE V) 



(7) 



In 2 



(ft)(m e c 2 f 



y F [Z, E v + Q n p] El (E v + Q np ) ((E v + Q np ) 2 - m e c 2 j [S v .][l- S e -]dE v , 



(8) 



In 2 



J F [Z. Q n p - E v \ El (Qnp - E v ) [{Qr, 



{ft)(m e c 2 f Jo 



E,, 



m e c 



1/2 



[l-Sv t ][l-S e -]dE v , 
(9) 



In 2 



(/t)(m e c 2 ) Jo 



F [Z, Q„p - £„] £ 2 (Q np - £„) ((Q np - E v ) 2 



1/2 



[SpJ^-Jd^, (10) 



r 



where E e and -Ej, are the appropriate electron/positron 
and neutrino/antineutrino energies. In these expressions 



the neutron-proton mass difference is Q 7 



1.293 MeV. 



Here ln2/(/i) is proportional to the effective weak cou- 
pling applying to free nucleons with (ft) the effective 
ft- value defined in Ref . [4y] . The weak matrix element is 
In 2/ (ft) oc G^r(l + 3g\), where Gf is the Fermi constant 
and gA is the ratio of axial to vector coupling for the free 
nucleons. In the BBN calculation the value for ln2/(/f) 
is normalized by the free neutron decay lifetime at zero- 
temperature. Here F [Z, E e ] is the relativistic coulomb 
correction factor (or Fermi factor) (46j . 



F(±Z,w) « 2(1 + s)(2pR) 



2(s- 



T(s + iff) 



T(2s + 1) 



(11) 



In this expression the upper signs are for electron emis- 
sion and capture, the lower signs are for positron emis- 
sion and capture, s = [1 — (aZ) 2 ] 1 / 2 , Z is the appro- 
priate nuclear charge (which is Z = 1 for the proton), 
a is the fine structure constant, r\ = ±Zw/p, and R 
is the nuclear radius in electron Compton wavelengths. 
R » 2.908 xl0" 3 A 1 / 3 -2.437A- 1 / 3 where A is the nuclear 
mass number and uj = (p 2 +m 2 ) 1 / 2 with m e the electron 



rest mass. This expression appears in the phase space in- 
tegrand of the weak rates which require a Coulomb factor 
in either the initial or final state 1">. lo If. 



7+ 



and S, 



are the phase space occupa- 
tion probabilities for electrons/positrons and neutri- 
nos/antineutrinos, respectively. For example, the 
[1 — S V J\ factor in X e - p is the Pauli phase space block- 
ing factor for processes which create a neutrino. In the 
limit that the neutrinos have thermally- shaped Fermi- 
Dirac distribution functions, these phase space occupa- 
tion probabilities become two parameter functions: 



S B . = 



1 



1 



(12) 



(13) 



The two parameters, T v and r\ Vt , correspond to neu- 
trino temperature and degeneracy parameter (the ratio 
of chemical potential to temperature), respectively. For 
example, a thermally-shaped neutrino phase space occu- 
pation probability function is graphed in Fig. Q] as the 
upper black curve. 



4 



The total weak neutron destruction rate is X n = X Uen + 
X e +n + ^n-dccay and the corresponding total weak pro- 



ton destruction rate is X p 
convenient to define 



Xj? e p + X e -„ + Xj) 



Atot — A n + X p . 



It 



(14) 



With this definition, the rate of change of the net elec- 
tron number per baryon, Y e , with Friedmann-Lemaitre- 
Robertson- Walker (FLRW) time-like coordinate t in the 
early universe will be 



dt 



= X n - Y P A t . 



(15) 



At early times where temperatures are high, the for- 
ward and reverse rates of these reactions are fast com- 
pared to the expansion rate of the universe. In this 
regime the neutron-to-proton ratio is just 



n 
P 



Xv c p + A e 



Xi/ e n T X e + n T X n decay 

This can be approximated as 

Ap c rj + A P 



n 
P 



•g p 



Xu e n T X e + n 



(16) 



(17) 



because neutron decay and the reverse three-body reac- 
tion are negligible by comparison at high temperatures. 
When the neutrino distribution functions have thermally- 
shaped Fermi-Dirac forms, the neutron-to-proton ratio is 
given by 



n ^ (A e - p /A e+n )+e-^+"'-g 

P \X e -n/X e + 7 



(18) 



where r\ Ve = fj, Ve /T is the electron neutrino degeneracy 
parameter, n e = (X e /T is the electron degeneracy param- 
eter, and £ is the neutron-proton mass difference divided 
by temperature, £ = (m n — m p )/T\j%. This equation 
is generally true whenever the lepton distribution func- 
tions have Fermi-Dirac forms and identical temperature 
parameters and whenever we can neglect neutron decay 
and its reverse process. Of course, at lower temperatures 
the neutrino and electron-photon plasma temperatures 
will differ and free neutron decay will be important. 

If the weak reactions occur rapidly enough to maintain 
chemical equilibrium, then the Saha equation, fi Ve +fi n = 
fi e - +/i p , can be used to predict the neutron-to-proton ra- 
tio. Interestingly, both the Saha equation and the steady 
state rate equilibrium condition in Eq. (|18[) , with the full 
lepton capture rates of Eqs. (fS ffTO)) . can be written as[5j 



(^ e -/ij/ e -5m np )/T 



P 



(19) 



This equilibrium neutron-to-proton ratio is shown in 
Fig. [5] as the dashed (green) line for zero electron and 
neutron chemical potentials, u e = /z„ e = 0. 
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FIG. 2: The neutron to proton ratio, n/p, as a function of 
temperature for three nucleosynthesis scenarios. The lower 
solid curve is for BBN with degenerate neutrinos and no neu- 
trino transformation, where L Vc — L Vt — L v = .05. The 
upper solid curve is the n/p ratio with the same lepton num- 
bers as above but now including a particular active-sterile 
neutrino transformation scenario. The dotted cure is the n/p 
ratio for standard BBN (no lepton numbers or neutrino oscil- 
lation). The dashed line is the n/p equilibrium prediction for 
standard BBN (no lepton numbers or sterile neutrinos) with 
enforced weak chemical equilibrium. 



As the universe cools, the weak reaction rates become 
slow compared to the expansion of the universe and the 
neutron-to-proton ratio falls out of equilibrium. This is 
called "weak freeze-out" and occurs over a range of tem- 
peratures. Fig. [5] shows the actual neutron-to-proton ra- 
tio evolving as a function of temperature for the stan- 
dard BBN scenario (thermal neutrino distribution func- 
tions and zero chemical potentials fj, e = [i Vil =0). At 
high temperatures, the actual neutron-to- proton ratio 
follows the equilibrium value and then around 1 MeV, 
the weak freeze-out commences. This happens because 
the weak rates have a stronger dependence on tempera- 
ture than does the expansion rate of the universe. The 
lepton capture/decay rates given in Eqs. ([5T fTU)l scale very 
roughly as T 5 (see Ref.[49] for the detailed temperature 
dependence) , while the expansion rate of the universe is 
ccT 2 . As a result, the neutron-proton weak interconver- 
sion rates eventually will fall below the expansion rate. 

Although the weak rates become relatively slow, they 
still have a significant effect on the neutron-to-proton 
ratio, even for temperatures well below T = 0.8 MeV. In 
fact, free neutron decay continues to lower the n/p ratio 
until there are virtually no more free neutrons or until the 
neutrons are sequestered in alpha particles, where they 
are effectively shielded from the weak interaction. This is 
illustrated in Fig.[2]where the dotted (blue) line continues 
to decrease until T « .08 MeV (when the neutrons have 
been captured during rapid alpha particle formation). It 
is important to correctly calculate the weak reactions in 
order to appropriately track the n/p ratio. This ratio 
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sets the scale, in varying degrees, for all the primordial 
element abundances [1J, k4| . 



III. NEW BBN CODE 

A nucleosynthesis code was written by Robert V. Wag- 
oner in 1969[52, 53] to track and time evolve the nuclear 
abundances and the neutron-to-proton ratio in an ex- 
panding cooling universe. It was later updated and re- 
vised by Lawrence Kawano in 1988 

This code time-evolves three main quantities, the elec- 
tron fraction, Y e , the baryon-to-photon ratio, 77, and the 
temperature, along with the primordial element abun- 
dances. It follows 48 nuclides using a reaction network 
composed of 168 nuclear reactions, whose rates have pri- 
marily been based on, and in some cases extrapolated 
from, laboratory cross sections. The main numerical 
technique is a 2nd order Runga-Kutta routine. 



The code also tracks the neutron-to-proton ratio by 
calculating the weak reaction rates using the stan- 
dard thermally-shaped Fermi-Dirac neutrino distribution 
functions, setting S Vc and So e as given in Eq. (|12[) and 
Eq. (HU). 

In their approach, electron energy is used as the inte- 
gration variable, instead of neutrino energy as given in 
Eqs. ([5l fTU|l above. To save computational time, they cal- 
culate only the sum of each of the forward n — > p rates 
and the reverse p — > n rates: 

A n = ^u e +n^p+e- + Ki+e+^p+P c + \i^>p+e~ +P„ (20) 
Ap Ap4_ e - -^u^+n "I" ^v e +p—*n+e+ ^p+e~+B K —*n' (21) 

With an algebraic trick, this simplifies the calculation by 
condensing the six phase space integrals (for each weak 
reaction rate) into two integrals: 



In 2 

(/t)(m e c 2 ) E 



x / EAE, 

J m e c 2 



(El - (m e c 2 ) 2 ) 



1/2 



(E e + Qnp) 



\E e Qnp) 



( e E e /T + I) t e -(E e +Q np )/T v -r,„ e + ]_) ( e -E e /T + ^ f e (E e -Q np )/T v - Vv 



(22) 



dE„. 



In 2 

(ft){m e c*f 



x J ^ 



(El~{ m ^) 2 ) 



1/2 



(E e + Q r 



( e E e /T + 1} ( e (B c +Q„ P )/T„+r,„ e + ]\ (eEe/T + ^ f e (Q np -E e )/T u + v „ e + ^ 



(23) 



dE P 



This algebraic trick requires the approximation of 
thermally-shaped Fermi-Dirac neutrino and antineutrino 
distribution functions. This summed rate cannot prop- 
erly treat the Coulomb correction, F[Z, E e ] , which should 
be included in the phase space integral of reaction rates 
which have an electron and proton in either the final or 
initial state. 

We have modified the Kawano/ Wagoner BBN code so 
that it can accommodate and integrate any arbitrary neu- 
trino and/or antineutrino distribution function with any 
specified time dependence. The majority of our changes 
lie in the weak reaction rate calculation. 

We first separated the summed neutron destruction 
and production rates, X n and X p . This enabled us to 
use non-thermal distribution functions and to change the 
neutrino and antineutrino distribution functions inde- 
pendently. Then, we removed a series approximation for 
A n and X p which is applied when the lepton numbers are 
zero. This approximation results in an erroneous ~ 0.5% 
increase in the neutron-to-proton ratio [Hil [55j. Further- 



more, we added the capability to separate a weak rate 
calculation into an arbitrary number of neutrino energy 
bins. This is useful for calculating a reaction rate where 
the neutrino energy spectrum is comprised of different 
functions over different energy ranges. 

For example, in Fig. [3J we have shown two electron 
neutrino distribution functions. The upper curve is just 
the standard thermally-shaped Fermi-Dirac distribution 
function, 

1 E 2 

fM) = (rK) eEv lT Va \„ a + r ( 24 ) 

which is consistent with the occupation probability de- 
rived from Eq. (|12|) . The lower curve is a distribution 
function resulting from a particular active-sterile neu- 
trino oscillation scheme described in Refs. [f| [HI]. In this 
scheme, electron neutrinos have been completely con- 
verted into steriles at low and high energies (1 and 3), 
leaving only active neutrinos in the center (2) energy 
band. To calculate a rate using this non-thermal dis- 
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FIG. 3: Two example electron neutrino distribution functions, 
where the upper black line is the standard thermal spectrum 
and the lower red line is a spectrum resulting from a partic- 
ular scenario for active-sterile neutrino mixing. The vertical 
dashed lines show where a weak rate calculation employing 
the lower distribution function would be broken up to be in- 
tegrated piece-wise in our new version of the code. 
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FIG. 4: Flow chart for our modified BBN calculation. 



tribution function, we break up the rate into three parts. 
The first part integrates from zero to t\ using the neu- 
trino distribution function f{E v /T) = 0. The second 
part integrates from t\ to ti using the modified function 
shown in 2. The third part integrates from 62 to 00 and 
again use f{E v /T) = 0. Finally, the total rate is calcu- 
lated by summing all three pieces. 

To perform these non-thermal piece-wise calculations 
in the BBN code, we completely replaced the original 
weak rate calculation with a series of four modules. These 
modules allow the user to define the distribution func- 
tions, break up the integration into specifiable pieces and 
define the energy ranges for each piece, and set any de- 
sired time/temperature dependence of the distribution 



functions. A flow chart of the weak rate calculations is 
shown in Fig. 0J At each time step, the BBN code calls 
the weak rate calculation subroutine, Module 1 in Fig. |4j 
to time-evolve the neutron to proton ratio and, subse- 
quently, all the nuclear abundances. 

Module 1 acts as the central line of communication in 
that it calls the other modules and reports back the value 
of the weak rates at every time step in the BBN code. 
In this module, the user can first define how many pieces 
to split the rate integration into for reactions involving 
either neutrinos or antineutrinos or both. For example, 
if the user wanted to use the lower non-thermal neutrino 
distribution function in Fig.[3]and a thermal antineutrino 
distribution function, the user can specify that the rate 
integrations involving neutrinos should be integrated in 
three parts and that rates involving antineutrinos should 
be integrated with one energy bin. 

Next, Module 1 calls Module 2 to retrieve the inte- 
gration limits for each piece, i.e., where the user wants 
each energy bin to begin and end. In Module 2, the user 
can define these integration limits and couple them to 
any time dependences desired. Module 1 makes an array 
with these limits so they can be accessed later in the inte- 
gration. This procedure can be extended to an arbitrary 
number of energy bins for any neutrino type. 

The first module calculates all six weak reaction rates 
by utilizing two main loops. These loop over the num- 
ber of energy bins. One loop calculates the two reaction 
rates that include neutrinos and the other loop calculates 
the four remaining weak reaction rates that include an- 
tineutrinos. The number of iterations for each loop is 
determined by the number of energy bins. Each loop it- 
eration integrates the weak reaction rates over the range 
of energy and neutrino distribution function specified for 
that energy bin. At the end of the iteration, each rate is 
summed. 

For every loop cycle, the first module calls the inte- 
grator which inputs the function to be integrated and 
the limits of the energy bins (from Module 2). The ma- 
trix elements and integrands for the six weak reaction 
rates, as shown in Eqs. ([ MTU)) , are retrieved from Mod- 



ule 3. Here, the electron occupation probability is set as 
S e = l/(e E °/ T + 1) and the neutrino and antineutrino 
occupation probabilities are called from Module 4. 

The sole purpose of Module 4 is to house the neutrino 
and antineutrino occupation probabilities. This makes 
it easy for a user to modify the neutrino distribution 
functions - by inputting analytic functions for S v ^ and 
<Sp e ~~ without having to modify any other portion of the 
weak rate calculation. The user can also define differ- 
ent functions or populations for each integration energy 
bin. After each energy bin is integrated, the total rate 
is summed and the values for the six weak reaction rates 
are returned to the main BBN code driver. 

Our modified Kawano/ Wagoner BBN code - which can 
now accommodate and integrate any arbitrary neutrino 
and/or antineutrino distribution function with any speci- 
fied time dependence - will be available to the community 
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T (MeV) 

FIG. 5: The rate of electron neutrino capture on a neutron 
as a function of temperature. The upper curve is Ai/ E n in the 
lepton number only case for lepton numbers of L„ e = L Vt — 
L Ufl = .05. The lower curve is the rate when there is active- 
sterile neutrino transformation along with the same lepton 
numbers as above. 

at bigbangonline.org[56j]. 



IV. EXAMPLE CODE RESULTS 

We have utilized this code to study nucleosynthesis 
abundance yields in the presence of a light-mass sterile 
neutrino over a range of lepton numbers [a HH • The lower 
red line in Fig.[T]shows a final non-thermal neutrino occu- 
pation probability function that can result from active- 
sterile neutrino transformation. In this particular sce- 
nario, we started with normal thermal electron neutrino 
and antineutrino distribution functions and an assumed 
initial lepton number. The lepton numbers that we have 
taken are within the range which is allowed by conven- 
tional BBN (primordial 4 He) considerations. But, of 
course, the point is that a sterile neutrino which mixes 
with an active neutrino can result in non-thermal neu- 
trino and/or antineutrino energy spectra which produce 
BBN abundance yields which can be quite different than 
in the standard scenario. This, in turn, could provide 
new, more appropriate constraints on lepton numbers 
or on active-sterile neutrino mass and mixing parameter 
space or on both. 

The presence of a significant net lepton number can de- 
lay significant sterile neutrino production until after the 
weak decoupling temperature. With a positive net lepton 
number, a Mikheyev-Smirnov-Wolfenstein (MSW) reso- 
nance occurs first for low neutrino energies. This reso- 
nance subsequently sweeps to higher neutrino energies as 
the universe expands and cools. At first, this resonance 
sweep process occurs adiabatically, efficiently converting 
all active neutrinos into sterile neutrinos. This contin- 
ues until the rate of active-sterile conversion becomes too 



fast to maintain adiabaticity. At this point, production 
becomes inefficient. However, at high enough resonance 
energies transformations can occur adiabatically again. 

Accurately following such a scenario requires all the 
modifications in our new code. Without being able to in- 
clude a dynamically changing neutrino distribution func- 
tion, for example, we could not calculate correctly the 
neutron-to-proton inter-conversion rates. In fact, in the 
example scenario presented here, not only are there non- 
thermal neutrino distribution functions to handle, but 
these change on time scales which are important to BBN. 
In Fig. we show the rate for electron neutrino capture 
on a neutron, the forward process in Eq. [21 as a function 
of temperature. The top curve is the rate when there 
is no active-sterile neutrino oscillation. The lower curve 
shows the decreased rate when there is active-sterile mix- 
ing and the final neutrino distribution function is that 
of Fig. [TJ By reducing the number of electron neutri- 
nos available for capture on neutrons, the capture rate 
is decreased. Additionally, the altered neutrino distribu- 
tion function also results in a modestly increased reverse 
rate (electron capture on protons). The depleted elec- 
tron neutrino distribution function in this scenario has 
the effect of increasing the electron capture rate because 
of the smaller neutrino phase space blocking factor. 

The final integrated effect in this scenario can be 
gauged by the changes in the light element abundances. 
For example, with adopted lepton numbers of L Ve = 
L v = L Ur = 0.05, which corresponds to a electron, 
mu, and tau neutrino degeneracy parameters of, rj Ue = 
r] v = t] Vt w 0.073 (i.e., near the conventional BBN up- 
per limits on these quantities), we see a 4.9% increase of 
4 He over the standard (no neutrino mixing and no lep- 
ton numbers) BBN value and a 12.7% increase over the 
4 He calculation with only lepton numbers included but 
no active-sterile neutrino oscillation effects. With this 
example scenario we find an increase in D/H (deuterium 
abundance relative to hydrogen) of 2.8% over the stan- 
dard BBN calculation and an increase of 6.9% from the 
lepton number only calculation. 

The increase in helium for these adopted parameters is 
likely unacceptable, exceeding observational bounds (57l 
|58| . [59| . Likewise, if the observationally-determined value 
of D /H can be increased in precision sufficiently (to bet- 
ter than ±5% .15]), it may be possible that D/H could 
compete with helium as an avenue for constraint of new 
neutrino physics. Ultimately, allowing for dynamically- 
altered neutrino and antineutrino distribution functions 
could add a new dimension to the way in which BBN and 
light element abundances might constrain new physics in 
the weak sector. 

We have also used our new code to apply a relativistic 
version of the Coulomb correction into the appropriate 
weak rate integrands]^ ■ This has never been done be- 
fore in the Wagoner/Kawano BBN code. 
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V. CONCLUSION 

We have developed an approach to Big Bang Nu- 
cleosynthesis (BBN) calculations where we can treat 
arbitrarily-specified energy distributions for all neutrino 
types, including v e and v e . We can also allow these distri- 
bution functions to be altered dynamically and follow all 
nuclear and weak reactions self-consistently with these 
alterations. This new approach can extend the useful- 
ness of BBN predictions for exploring and constraining 
new physics in the neutrino and weak interaction sectors. 

Examples of such new physics include active-sterile 
neutrino mixing and particle decays that have neutrinos 
in the final state. We have given an explicit example of 
the former scenario. In this example we have demon- 
strated how active-sterile neutrino oscillation physics 
can alter neutrino or antineutrino distribution functions 
on short time scales, alter the neutron-proton inter- 
conversions rates, and so modify BBN abundance yields 
over those of the standard scenario. 

Our calculations hold out the promise that light ele- 
ment abundances could place the best constraints on pri- 
mordial lepton numbers and active-sterile neutrino mix- 
ing parameters when the sterile neutrino mass is in the 
~ 1 cV range. Present laboratory experiments, like mini- 
BooNE, are sensitive to neutrino flavor mixing in the 
active-sterile channel at the ~ 1, eV mass scale only when 
the appropriate effective 2x2 vacuum mixing angle sat- 
isfies sin 2 29 1CD 4 . By contrast, in the presence of 
a net lepton number, BBN abundance yields might be 



significantly altered for active-sterile neutrino mixing pa- 
rameters for sin 2 29 > I CP 8 . The greater reach in vac- 
uum mixing angle afforded by BBN considerations stems 
from: (I) the long (gravitational) expansion time scale 
of the early universe which dictates the MSW resonance 
sweep rate and sets the minimum mixing angle required 
for adiabatic and efficient conversion of the active neutri- 
nos into sterile species; and (2) the significant sensitivity 
of the neutron-proton weak inter-conversion rates to al- 
terations of the neutrino or antineutrino energy distribu- 
tion functions. Our new calculations allow us to follow 
simultaneously and self-consistently both of these effects 
along with all relevant weak, electromagnetic, and strong 
nuclear reaction rates. 

This new approach is incorporated into an update of 
the Kawano/Wagoner BBN code - which can now ac- 
commodate and integrate any arbitrary neutrino and/or 
antineutrino distribution function with any specified time 
dependence. We will soon make this code available to the 
community at bigbangonline.org. 
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